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Abstract 

Direct numerical simulations (DNS) of passive scalar mixing in isotropic turbulence is 
used to study, analyze and, subsequently, model the role of small (subgrid) scales in the 
mixing process. In particular, we attempt to model the dissipation of the large scale (su- 
pergrid) scalar fluctuations caused by the subgrid scales by decomposing it into two parts: 
(i) the effect due to the interaction among the subgrid scales, and, (ii) the effect due 

to interaction between the supergrid and the subgrid scales, Model comparison with 

DNS data shows good agreement. This model is expected to be useful in the large eddy 
simulations of scalar mixing and reaction. 
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1 Introduction. 


Large eddy simulation (LES) of turbulent mixing and reaction is a possible means of cal- 
culating complex flows not amenable to the full direct numerical simulations (DNS) [1], [2]. 
In performing the LES of combustion, one is faced with a modeling problem not present in 
the LES of Navier-Stokes equation, that of molecular mixing of scalars. Many combustion 
problems of practical interest are of the non-premixed type. In non-premixed combustion, 
scalar mixing, to a large extent, controls the rate of reactant conversion. Scalar mixing is pre- 
dominantly a small-scale phenomenon occurring in scales not resolved in LES calculations. 
Adequate modeling of the subgrid-scale mixing process is, hence, crucial to the success of 
LES of reacting flows. 

Consider the mixing and reaction of a scalar field <j>(x, t) in an incompressible turbulent 
velocity field u(x,t): 


d<f> d<f> d 2 <f> 
dt +Ua dxl = + "' W ’ 


(1) 


where // 0 and w are the molecular diffusivity and the chemical conversion rate of the scalar. 
(Repeated indices imply summation.) In a turbulent flow, u a and <j> fields are made up of 
the entire gamut of scales ranging from the energy containing scales to the dissipative scales. 
In a high Reynolds number flow of complex geometry, a DNS calculation that resolves all of 
the length and time scales represents a formidable task. Large eddy simulation entails the 
calculation of the energy containing (called supergrid or resolved) scales only. The dissipative 
and dispersive action of the small (called subgrid or unresolved) scales on the large scales 
are accounted for by using suitable models. 

The resolved field of any quantity A is obtained from its full field A by employing a filter 
G of specified width A: 


/ OO 

A(x',f)G(x-x',A)dV. 

-OO 


(2) 
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The unresolved field, A', is given by 


A'(x, t) = A(x, t ) - A(x, f ). 


( 3 ) 


When the filter is applied to the scalar evolution equation (1) we get 


dp dp 

dt 0 dx a 


H 0 


cPp 

dx a dx a 


+ w(<j>) - 


djUgP - UgP) 

dx a 


( 4 ) 


The last two terms of the right hand side (RHS) of the above equation need closure modeling. 
In the absence of reaction, the effect of the unresolved scales on the resolved scales manifests 
only via the last term on the RHS which is usually modeled using a eddy viscosity type 
model: 

—7 _T 

^0,4* — t*T r\ i \^/ 

OX a 

where the eddy diffusivity coefficient fij can be obtained using the dynamic subgrid scale 
modeling strategy outlined in [3] or the multiple scale elimination scheme [4]. 

The perils of modeling the filtered reaction term as 




(6) 


are well known. Following the procedure used in [2], this term is expressed exactly as 

yV>= 00 

w(x,t)= (7) 

Jy= — oc 

where ip represents the probability space variable corresponding to the scalar <p and the 
function / is the ‘large eddy pdf’ (probability density function). If the function / is known, 
equation (7) can be used to calculate the filtered conversion rate. There are two possible 
ways of obtaining /. The computationally simpler option, in the same vein as the assumed- 
pdf approach to turbulent combustion modeling, is to presume the form of f(ip] x,t) knowing 
its first few moments at each location [1]. Most assumed pdf’s require the specification of 
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the first two moments: in this case, <p and <f > 2 . The filtered field <j) is taken to be known; (f> 2 
may be determined by solving its governing equation: 


d<t > 2 , d<t > 2 d 2 <f > 2 dM 2 - xr a <i > 2 ) — - 

■qT + U aR = o R + 2 w(<f>)<f> - 2n 0 - — 

ot dx a dx a dx a dx a v ^ dx a dx a 


d<j> d<f> 


( 8 ) 


One of the terms that need modeling in the above equation (8) is the last term on the (RHS) 
called the filtered scalar dissipation £ 4 . The second, and perhaps the more accurate, option 
of obtaining the function / is to solve its governing equation [2]. The most important term in 
the large-eddy pdf governing equation that needs closure modeling is the filtered conditional 
diffusion. Most of the current models for the filtered conditional diffusion (eg., least mean 
square estimation, mapping closure methodology) require that the filtered scalar dissipation 
be known. 


It is clear, then, that the role of subgrid scales during mixing, manifesting as the filtered 
scalar dissipation, is a key process that needs closure modeling. Before attempting its model- 
ing, it is important to understand the basic difference between statistical mean of a quantity 
denoted by angular brackets (e.g., (<p 2 )) and its filtered counterpart represented by a over- 
bar ( q !> 2 ). The difference is most easily understood in isotropic turbulence where (<^ 2 ) and 
are cons tants in space, but their filtered counterparts exhibit spatial dependence. 
In Fourier space, the spectrum of the statistical mean is a delta function at zero wavenum- 
ber; whereas, the filtered quantity typically has a broad-based spectrum with the highest 
amplitudes in the small wavenumbers (corresponding to large scales of motion). The lack 
of a model that can adequately simulate the correct spatial (or, correspondingly, spectral) 
behavior of the filtered quantity may nullify any advantage that the LES may have over the 
traditional moment method of calculating mixing and reaction. The filtered scalar dissipa- 
tion model currently used in literature [1] is an adaptation of the mean scalar dissipation 
model used in traditional moment methods: 

c - W 

£<t> C<f, ^ , (9) 
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where C t, is a constant of order unity and r is the time scale of turbulence. We will first 
demonstrate that this model is inadequate for capturing certain important features of filtered 
scalar dissipation observed in DNS data. We will then proceed to derive a more adequate 
model expression. 

As a first step, in this paper, only inert scalar mixing is considered. In practically all 
combustion calculations to date, the form of the mean scalar dissipation model used is the 
same for reacting as well as inert flows. Following suit, even though the current model is 
developed only for LES of inert scalar mixing, it is recommended for use in reacting cases 
also. In the discussions so far, G was permitted to be any positive- definite filter. The 
derivation below is based on a spectral cut-off filter. We speculate that the functional form 
may be valid for other filters also since the objective of most filters is indeed to separate the 
large and small scales, even if only approximately. 

In Section 2, we perform DNS of passive scalar mixing. We study the DNS data to 
develop intuition into the behavior of ^ and The DNS data is also used to validate 

some modeling assumptions and ultimately to evaluate the model itself. Section 3 contains 
model development and validation. We conclude in Section 4 with a summary. 

2 Numerical simulations and scalar analysis. 

Direct numerical simulations of passive scalar mixing similar to that of Eswaran and Pope 
[5] are performed. The initially non-premixed scalar field evolves in isotropic turbulence 
that is held statistically stationary by small wavenumber forcing. The scalar field is initially 
permitted values of only 4-1 or —1 with equal probability. The calculations are performed 
on a 64 3 grid, and the Reynolds number based on the Taylor length scale is about 40. The 
results presented in the paper are for a scalar of Prandtl number unity. Simulations were 
also performed with non-unity Prandtl numbers and those results are consistent with the 
presented results. The evolution of the scalar energy {E#) and dissipation ( D $ = k E,p(k)) 
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spectra are presented in figures la and lb. respectively. As can be seen from the figures, 
due to the rather modest Reynolds number of the simulation, there is no inertial range or 
separation of scales. The total number of energetic wavenumbers in a 64 3 simulation with 
spectral-cut-off de-aliasing is only about 30. 

We are primarily interested in modeling the filtered version of the equation 

df_ df^_ (Pf_ d<fr 

dt a dx a ^ °dx a dx a dx a dx a 

Due to the presence of the gradient-squared term (last term on the RHS), the above equation 
is different from that of a passive scalar. A complete understanding of the spectral behavior 
of 4 i 2 and the gradient-squared term is important for modeling their filtered counterparts. By 
writing equation (10) in spectral space and multiplying it through by the complex conjugate, 
<f) 2 *( k), we get the spectral budget for d> 2 -energy: 

[j t + 2M 2 ]iv(k) = 27>(k) - 4S(k), (11) 

where, 

EM k) = ^(k)^(k^ (12) 

S( k) = HoReal[p-^-(k)<j> 2 *(k)] = Real^ k)<£ 2 *(k)]. 

OX q C/X a 

In the above equations, superscript * denotes complex conjugate of the quantity and ~ 
indicates Fourier transform. For the sake of convenience, we do not use ~ over d>(k), £${ k) or 
</> 2 (k). The second term on the left hand side (LHS) of the above equation (11), D^{ k) = 
2 fi 0 k 2 E 4 > 2 (k ) 1 represents the loss of <?i> 2 -energy due to molecular action. The first term on the 
RHS, 7>(k), represents the net spectral transfer of <^ 2 -energy into wavenumber k due to 
the cascading effect of the velocity field. The term S(k) represents <^ 2 -energy transfer into 
wavenumber k due to the gradient-squared term. 
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2.1 Data analysis. 

Behavior of E & ( k ). The isotropic one-dimensional spectrum E^(k) calculated from DNS 
data is given in figure 2a. Most of the <f> 2 energy is initially present at the large scales or 
small wavenumbers. With passage of time, larger wavenumbers become more energetic at 
the expense of the smaller wavenumbers. This transfer of <f ) 2 energy from large to small 
scales could be due to the effect of the velocity field cascade term (!>(*)) and/or the 
scalar gradient-squared term (S{k)). At the final stages of mixing, as is to be expected, the 
spectrum decays over the entire wavenumber range. 

Behavior of S(k). In physical space, the gradient-squared term always decreases the level 
of <j> 2 (equation 10) lowering its spatial mean. This means that, referring to equation (11), the 
value of S { k = 0) has to be positive. The behavior of S{k) at other wavenumbers is not clear. 
By depleting <f >‘ 2 at different rates at different physical locations, the gradient-squared term 
may be capable of creating or augmenting spatial fluctuations in the (j) field. If the gradient- 
squared term indeed creates or augments <f> 2 fluctuations at a wavenumber, then it would be 
a -energy source at that wavenumber. This means that, again referring to equation (11), 
the value of S(k) at that wavenumber would be negative. In figure 2b, the function S(k) 
(obtained by summing over all wavenumbers of a given magnitude k) is plotted as a function 
of k(> 0) at various stages of mixing. Initially, S(k) is negative everywhere indicating that 
the gradient-squared term acts as a ^-energy source at all wavenumbers k > 0, With time, 
the magnitude of this energy source decreases. At later times, S(k) is positive at small 
wavenumbers indicating that it acts as a energy sink at these wavenumbers while still acting 
as energy source at higher wavenumbers. At the very final stages of mixing, we speculate 
that S(k ) would would be negative everywhere, thereby reducing the fluctuations in O 2 at 
all scales. 
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DNS data vs. Moments-method model. The model for LES filtered scalar dissipa- 
tion given in equation (9), adapted from the moment-method model for the mean scalar 
dissipation, leads to the following behavior for S(k)\ 

S(k) = Realist 2 (k)cf> 2 *(k)\ = —E+*(k) > 0. (13) 

T T 

According to this model, S(k) is always positive, and, hence, a sink of </> 2 -energy at all 
wavenumbers and at all times. This, clearly, is inconsistent with the behavior of S(k) 
observed in the DNS. Even at high Reynolds number, there is no reason why the gradient- 
squared term cannot transfer 0 2 -energy from large to small scales and, hence, take negative 
values at large wavenumbers. Therefore, we conclude that the model for filtered scalar dissi- 
pation currently used in literature (equation 9) is inadequate for capturing some important 
features of scalar mixing. 


2.2 Decomposition of scalar dissipation. 

If we use the spectral cut-off filter, the scalar and the velocity fields can be decomposed into 
super and subgrid quantities as follows: 




and, similarly, 


i (]< t) = / if k < A c ; 

* [Kt) 10 \f k < A c . 


{^(M) i 


if k < A c 
if k > A c 


if k < A c 
if k < A c . 


Here, A c is the wavenumber that separates the resolved and unresolved scales and is typically 
located in the inertial sub-range. The turbulent velocity field u(k,f) is determined from the 
Navier-Stokes equation and the scalar field evolves according to equation (1) with no chemical 
reaction. 

We seek to model the filtered scalar dissipation 


£*(x) — Po 


d<f> d<j) 

dx a dx a 


(14) 
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for the case of spectral cut-off filter. The filtered scalar dissipation can be decomposed into 
three parts as follows: 

£<t> = £$ < + £$ y + £$ < 

d~4> d* ^ W l9 ~WW n -x 

^ °dx a dx a ^ °dx a dx a ^ °dx a dx a 

In the above equation (15), the first term on the RHS is the scalar dissipation due to 
interactions among the resolved scales. This term is closed in the supergrid variables and no 
modeling is required. The second term on the RHS of (15) represents the scalar dissipation 
caused by interactions among the unresolved wavenumbers. This term is always positive 
and requires closure modeling. The third term on the RHS represents the dissipation due 
to interactions between resolved and unresolved scales. This term, which also needs to be 
modeled, can be either positive or negative and is responsible for ‘backscatter or energy- 
transfer (of ^-energy) from small to large scales. If a sizable spectral gap exists between the 
resolved and unresolved wavenumbers, this term will be negligible. However, the turbulence 
spectrum is continuous on either side of the cut-off wavenumber. Analyses of DNS data have 
demonstrated that interactions of this type play a major role in spectral energy transfer (Zhou 
and Vahala [6]). Therefore, an accurate description of turbulent mixing must incorporate, 
as a fundamental building block, this interaction between the resolved and unresolved scales 
(Rose [7]). 

The roles of the three components of dissipation in the spectral budget of (j ) 2 -energy are 
now examined. In particular, the evolutions of S << (k), S >< {k ) and S >:> (k) - 

S«(k) = £<<(k)<f> 2 *(k) S><(k) = 2£><(k)^(k) S»(k) = £><{k)<p 2 *{k\ (16) 

the three components of S(k ) corresponding to the three components of dissipation - are 
investigated. The results from calculations using DNS data, with A c = 8, are presented in 
figure 3. At early times, 5 << (it) is much larger than the other two components. At later 
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times, the three components are of the same order of magnitude. In a typical high Reynolds 
number turbulent mixing problem of practical interest, we would expect that S << (k) would 
be much smaller than the other two components. Two other important observations, from 
figure 3, which we expect to be valid at high Reynolds number also are: (i) 5 >> is more 
important than S ><: at small supergrid wavenumbers; and, (ii) S >< is larger than 5’ >> at 
the largest supergrid wavenumbers. 

Our objective is to construct models for ££ K and ££> which will lead to the behavior 
of S ><: (k) and S >> (k) consistent with above observations for the supergrid wavenumbers 
(k < A c ). 

3 Model development. 

In a high Reynolds number scalar mixing problem with adequate scale separation, most of 
the scalar dissipation, £$( k), is due to £^ > (k) and it occurs at scales that are much smaller 
than the length scale corresponding to the cut-off wavenumber, A c . On the contrary, Rose [7] 
suggests, that most of the contribution towards the back-scatter term, £^ < (k < A c ), comes 
from the interactions between the supergrid and the largest subgrid scales, (wavenumbers 
whose magnitude are of the order A c ). In figure 4, we present the contributions of various 
sub-sets of the subgrid range wavenumbers towards S >< (k), calculated using DNS data. It 
is clear from the figure that most of the contribution indeed is from scales only slightly larger 
than the cut-off wavenumber. In fact, the contribution from wavenumbers larger than 2A C 
is negligible. These observations lead to the following modeling simplifications: only the 
near subgrid scales contribute towards ££<(&); and, the only contribution towards £?>(£) 
are from the far subgrid scales. The two components of dissipation are, then, modeled 
separately employing the assumption of a spectral gap in deriving the model for £$ > (k) 
only. To facilitate the implementation of this simplification in our model development, we 
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subdivide the subgrid wavenumbers into far and near subgrid scales: 


lt(k if A ; < < Ajv ; . HA c <k<A N 

(lc ’ t) “ to if fc > Ajv- ’ * (k ’'' U'(k,() if k > A/v. 

where, Aat (< 2A C ) is the cut-off between far and near subgrid scales. 

3.1 Model for £’^ > (k). 

In terms of these new variables we can write 


dd>' d<f>' d 4> y d <^ > 

£1 = ^o^—— ~ Mo- 


( 17 ) 


'dx a dx a dx a dx a 

In deriving the model for this term, we assume that the scalar field is composed only of 
the supergrid and far subgrid wavenumbers which are separated by a sizable spectral gap. 
Clarification of the term ‘spectral gap’ is called for here. By spectral gap we mean the 
gap between the supergrid (energy containing scales) and the scales which contribute most 
towards ^^(k). The scales of motion that contribute most towards this dissipation are of 
the order of the Kolmogorov length scale. Therefore, the larger the Reynolds number, the 
wider is the spectral gap in our definition, and hence, the more valid is the analysis presented 
below. 

It is known from direct interaction approximation [8] and renormalization group theory 
[9] estimates that, if there is spectral gap, then the effect of the small scales on large scales 
of motion can be reasonably well accounted for by using a gradient diffusion approximation. 
Therefore, the total scalar field evolution equation 


d<f> d<f> _ d 2 4> 

dt ° dx a dx a dx a ’ 


( 18 ) 


on elimination of the subgrid scales leads to the following evolution equation for the filtered 
<f> field: 

dt +Ua dx a dx a dx ( ) 
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In a high Reynolds number flow the turbulent diffusivity fi T is much larger than the molecular 
diffusivity fi 0 . In this paper, we are interested in modeling the evolution equation of the 
filtered (f > 2 field, whose total field evolves according to 


dt +Ua dx a 


d 2 ^ 2 d<j> d<f> 

^ 0 dx a dx a ^ 0 dx a dx a 


We propose that, much in the same way that equation (18) on filtering leads to equation 

(19), the total <b 2 evolution equation (20) would, on filtering, lead to 

dd > 2 _ d<f> 2 d 2 4> 2 d4> d(j> 

iH +u *dxi = _ 2i>T dr,dxi' ( 21 > 

Comparison of this modeled equation (21) with the exact evolution equation (8) (after setting 
w {4 > ) = 0) yields the following. The scalar flux model is 


d{u a (p 2 - u a <j> 2 ) ' N d 2 (j) 2 


The implied model for filtered dissipation is 


d(f> d<f> _ d<f) d(f> d ^ , 

^ °dx a dx a dx a dx a + dx a dx a 


dx a dx a 


dl 

dx a dx a 


This leads to 


c» - /- , d<f> d<f> d<f> d<j) 


The eddy diffusivity coefficient fij can be obtained using the dynamic subgrid scale modeling 
strategy outlined in [3] or the multiple scale elimination scheme [4], [7], 

As mentioned earlier in this subsection, this model for ££> may not be accurate un- 
less there is adequate separation between the supergrid (energy containing) scales and the 
dissipation scales (of the order of Kolmogorov scales). Hence, the above model cannot be 
expected to perform well for low Reynolds number mixing problems. With complete cog- 
nizance of this fact, we compare the model computations of S >:> (k) with DNS data in figures 
3, using 

Vt _ = 0 ) 

/Jo -£<<(k = oy (25) 
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In the model calculations, the function is of course taken from DNS data. Except 

at the final stages of mixing, the model performance is unexpectedly good. The model 
appears to capture the wavenumber dependence of the S >> (k ) reasonably well. The model 
performance can be expected to be better for higher Reynolds number mixing problems of 
practical interest. 


3.2 Model for ^^(k). 


Many of the shortcomings of LES calculations are generally attributed to the poor, or even 
lack of, model for the backscatter effect. Derivation of an adequate model for the backscatter 
dissipation represents a major thrust of this paper. Before we can model the backscatter 
effect of the near subgrid on the supergrid scales, the effect of the far subgrid on the near 
subgrid scales must be accounted for. We start with the assumption that the far subgrid 
affects the near subgrid in the same way that it affects the supergrid scales. In other words, 
the scope of equation (19) is now extended to include the near subgrid scales as well. Clearly, 
equation (19) is not as accurate for the near subgrid scales as it is for the supergrid scales, 
due to the lack of spectral gap between the near and far subgrid scales. Hence, in assuming 
that the near subgrid scales are also governed by the enhanced diffusivity equation, the 
backscatter effect between the near and far subgrid is ignored. We argue that, while this 
neglection may be of concern if we are attempting to model the near subgrid scales, it may 
be unimportant in the modeling of the interaction between the supergrid and near subgrid 
scales. Subject to the above simplification, the spectral evolution of the near subgrid scalar 
field is given by 


= —ik a j <i 3 ;u a (k - 


(26) 


where, i = All of the effects of the far subgrid scales are accounted for by enhancing 

the diffusivity from no to fij- From this point, the velocity and the scalar fields are considered 
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as being composed only of the supergrid and near subgrid wavenumbers. Subject to these 
assumptions we have 


= 2// 0 


d<f> d<j)' 
dx a dx a 


In spectral space, can be written as 


2/j.q 


d(j) d<f> + 

dx a dx a 


(27) 


^^(k) = -2h t (J d z j j a (k - j) a <f>( k-j)^+(j)). (28) 

At this stage, one has the option of either deriving a modeled evolution equation for 
£^1 or see k a n algebraic model. For the sake of simplicity, we will restrict ourselves to 
an algebraic expression for ££ < closed in terms of the supergrid quantities. Most algebraic 
subgrid scale models - e.g., the Smagorinsky model or the dynamic subgrid scale model for 
the velocity field — make the implicit assumption that the subgrid scales are in some kind 
of an equilibrium with the supergrid field, so that no evolution equation is needed for the 
subgrid scales. We make this equilibrium assumption explicitly and employ its implications 
to derive the model for £><. 

We start with the assumption that the scalar spectrum of the subgrid scales is in equi- 
librium with the large scales: 

dEJk) 

— — = 2T^(k) - 2fi T k 2 (f>(k)<f>*(k) « 0, for A c > k > A N . (29) 

In the above, summation is performed over all vector wavenumbers of a given magnitude. 
(Again following Rose [7], we treat the near subgrid scales as nearly isotropic, even if the 
supergrid scales are not.) The scalar energy transfer into a wavenumber k due triadic 
interactions between the velocity and the scalar fields, denoted by T^k), has the following 
form: 

7V(k) = Real[-ik a J d 3 ju a ( k - j)<£(j) <£ + *(k)]. (30) 

The scalar energy dissipated in wavenumber k is Z^( k) = -2n T k 2 <f>(k)<f>*(k). The simplifying 
assumption in equation (29) is now verified using DNS data. In figure 5, we plot T^k) and 
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D<t,(k) for the subgrid wavenumbers (k > A c = 12) at various stages of scalar field decay. 
At initial times, T^k) is slightly larger than D+(k) indicating the build up of energy in the 
near subgrid scales. At later times, due to the decay of the spectra at all wavenumbers, 
D^k) is slightly larger than T^k). However, throughout the mixing process, the difference 
between D^k) and T+{k) is much smaller than their absolute magnitude confirming that 
the simplification in equation (30) is quite reasonable. 

Based on the observation that D 4 (k) « T*(fc), we suggest that, in an average sense, there 

is balance at each wavenumber vector: that is, 

D*(k) « r*(k) ( 31 ) 

^ T lb 2 [^ + (k)]0 + ’‘(k) « Real[-ik a j d 3 ju a (k - j)<Mj)<£ + *(k)]. 

If the above assertion is true, then we must have 

<j> + {k 1 t) = ^ 3 i u «( k _ + c 'o^ + (k,t)], (32) 

where, Co is a constant and <ji + (k) is such that 

<£ + (M)<£ + ’(M) = 0. ( 33 ) 

Substituting equation (32) into the expression for (k) in equation (28) we get 

£><(k) = -2pt( j d 3 j Ja(k ~ - })-^{~ijb J d 3 j'u b ( j - j'W) + Co^ + (j)p4) 

= 2 (J J d 3 jd 3 j'C-p)j a (k -j) a <f>(k - j)«&(j 

-2Co( j d 3 jj a {k - j)a(4)^( k - j)^ + (j)>- 

In the above equation, the integration over j is restricted to subgrid wavenumbers, whereas, 
the integration over j' is over all wavenumbers. Considering that the subgrid now consists 
only of a narrow band of scales, we propose that, to leading order, 

( j J d 3 jd 3 j'j a {k ~ j)aCjj) ^(k - j)u 6 (j - jW)> ( 35 ) 

~ (J J d 3 jd 3 j'j a (k — j) a (-^)4>(k — j)«fc(j — j j ))• 
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After this truncation to the leading order term, the only remaining unclosed term in the 
model for ££ < is the one containing 4> + . Calculations of this term from DNS data (result 
not shown) appear to indicate that it is not small in magnitude compared to the other term 
(equation 32) in the model expression for £>< . However, we know from equation (33) that 
this term does not affect the level of the scalar energy. Since our main interest in £% < is to 
be able to calculate the level of scalar energy of the supergrid (<^> 2 ) correctly, we suggest that 
the term containing d> + be dropped from the model expression. In other words, we suggest 
that exclusion of this unclosed term would still lead to the correct behavior of the ultimately 
important quantity S >< {k) (see equation 16 for definition), even if the behavior of £><(&) 
itself may not be quite correct. The extent of validity of this exclusion will be borne out 
when the model calculations of S y< (k) are compared against DNS data. 

The expression for ^k) (equation 34) can now be approximated as (the angular brack- 
ets indicating averaging are dropped because only the supergrid quantities survive our sim- 
plification) 

£><(k)*2j (36) 

The above equation represents a closure model for the backscatter dissipation in spectral 
space. Most of the LES computations are performed in physical space, and, therefore, a 
closure model in physical space would be very useful. Inverse Fourier transformation of 
equation (36) would lead to an expression that involves Greens integral over space, and, 
hence, not local in space. The non-locality of the closure model for ££< will preclude it from 
being useful in LES computations of high Reynolds number flows of complex geometry. 

A model local that is local in physical space can be obtained if we assume that the contri- 
bution towards £^ < comes only from a thin shell of subgrid wavenumbers, an approximation 
that is justified by the finding in figure 4. Then, in equation (36) we can set |j| « A c (recall 
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that the dummy wavenumber j is restricted to near subgrid scales), leading to 


^^(k) « ^ J d 3 jd 3 j'(k - j) a ja{ijb)ub{j - ? ,t)<j>(k - j ,*)• (37) 

By performing a convolution integral over j' (which is over the entire wavenumber domain), 
we get 


(k) = ^ J d 3 j(k - j)J„(ij : i)M](j,()^(k-j,() 

= _1 U °j— , 

A 2 J J dx a dx b 


(38) 


d 2 u b <f) d<j> 

(k-j.O- 


The integration over j (subgrid scales only) can be split into two parts 


J subgrid dx a 8x b 8x a Jtotal dx a 8x b 3x a 


■ 8 2 u b (j> ^d(f> 


(39) 


-L 


d 3 J 


d 2 u b <f> d<j> 


supergrid 8x a () X b 8x a 




Each integration can now be separately evaluated: 

f * • d 2 W r .x /. ; .x _ r 

Jtotai d J dx a dx b ^dx a ’ dx a dx b dx a 


](k) 


(40) 


f ,3 -Pubdr^di . ,, f i\H!L{ir- i /I = i_22Sf-.22.inc 

S, m J 1 d^y ) wj ' ] k ' J ’ <) “ L,.i ij dx,dx b b’ ) dx} k J,i) [ at.ax»*r. 1( ’ ’• 


. d 2 u b (f> 


dj> 


d 2 u b <j> d(j> 


Therefore, the exact inverse transform of equation (38) yields, 


^ 2 dcf>.d 2 u b (t) d 2 u b 4> 

[X) ~~A 2 dx a [ dx a dx b dx a dx b l 


(41) 


The parameter A c can be expressed in terms of the mechanical dissipation rate (£) and the 
kinetic energy (K) of the subgrid scales (e.g., Rubinstein and Barton [10], Zhou et al. [11]): 


Ac 


£ 

K 3 / 2 


(3C*/2) 3/2 , 


(42) 


where Ck is the Kolmogorov constant. To obtain the relations in equation (42), the velocity 
field is assumed to be stationary and the small scale autocorrelation is taken to be determined 


16 



by the Kolmogorov spectrum. After substituting A c from equation (42), the model expression 
for £ ^ in terms of the subgrid turbulent velocity field variables is given by 


£><( x ) - r>< — f 

* * £ 2 dx a dx a dxb dx a dxb 


(43) 


where is independent of the wavenumber. 

We now compare our model against DNS data. We reiterate that the main emphasis of 
the modeling effort is to be able to calculate the correct spectral behavior of the supergrid 
component of 0 2 which evolves according to equation (8). To examine if the proposed 
model accomplishes the stated objective, we compare the model calculation with DNS for 
the quantity S >< (k ) = ££ < <f> 2 *(k) in figure 6. The comparison is performed at four different 
stages of the mixing process with A c = 12; C>< is taken to be that numerical value which 
gives the best agreement at each time. The model agrees very well with DNS data, especially 
at early times. As required, the model appears to capture the wavenumber dependence 
of S(k) very well. The modeling is complete except for the specification of C>< which 
determines only the overall level of S(k) and not its spectral character. The determination 
of is deferred to a future effort. 


4 Conclusion. 


In this paper, we first study the behavior of filtered scalar dissipation - arising in LES 
modeling of turbulent scalar mixing - with the help of DNS data. Then, we employ the 
insight gained to develop a model expression for the filtered scalar dissipation. It is found 
that its accurate modeling requires that two physical phenomena be accounted for: (i) 
the dissipation due to interaction among the subgrid scales, denoted by £?>; and, (ii) the 
dissipation due to the interaction between the resolved and unresolved scales, given by £><. 

Arguing that the dominant contribution towards £% > comes from the interaction among 
far-field subgrid wavenumbers, we introduce an artificial spectral gap in the derivation of the 
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model for this term. The model for ££> turns out to be an enhanced molecular diffusion 
type term, as it should be in the distant interaction approximation limit. Although the 
assumptions leading to this model are valid only for high Reynolds number turbulent mixing, 
comparison with low Reynolds number DNS data shows surprisingly adequate performance 
by the model. 

It is also demonstrated using DNS data that most contribution to the ‘backscatter- 
dissipation, comes from a narrow band of subgrid wavenumbers only slightly larger 

than the largest supergrid wavenumber. After establishing that the scalar spectrum of this 
near subgrid wavenumber band is in quasi-equilibrium with the large (supergrid) scales, we 
derive the model for in its most accurate form in the spectral space. Inverse transform 
into physical space leads to a non-local model, which in large scale LES calculations, could 
lead to excessive computational requirements. A less accurate, but computationally more 
feasible, model for this term that is local in physical space is derived by assuming that only 
a thin shell of wavenumbers contribute towards the ‘backscatter’. Comparison of this model 
with DNS data shows good agreement. 

The present model appears to be a clear improvement over the previously used model 
for filtered scalar dissipation. The previous model, which is a literal extension of the model 
for the statistical average of scalar dissipation used in moment methods, makes no use of 
the additional information about the supergrid scalar field available in a LES computation. 
The present model makes use of the additional supergrid scalar field information available, 
leading to a more accurate model. Owing to the reasonably sound theoretical foundation 
of the new model, coupled with its demonstrated ability to capture important features of 
DNS mixing data, we expect that the model derived in this paper would be very useful for 
accurate LES computations of scalar mixing and reaction. 
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Fig. la: Scalar Spectrum 



Wavenumber k 


Fig. lb: Dissipation Spectrum 



Wavenumber k 


Figure 1: Evolution of scalar spectra calculated from DNS data: (a) Energy spectrum; 
(b) Dissipation spectrum. Times: tl = 1.25 r,,; t2 = 2.5 r,; t3 = 3.5 r and t4 = 5.25 
where the Kolmogprov time scale of turbulence ~ 0.1. 
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Fig. 3c: Time t3 



Fig. 3d: Time t4 



Figure 3: Evolution of the various components of S(k) calculated from DNS data. Evo- 
lution of S >:> {k) as computed by the enhanced diffusion model is also presented. Cut-off 
wavenumber A c = 8. Times same as given in figure 1. 
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Fig. 5a: Time tl 





D^k), T^k) D t (k), T t (k) 


Fig. 5c: Time t3 



Fig. 5d: Time t4 



Wavenumber k 


Figure 5: Evolution of D(k) and T(k) of the scalar spectrum evolution equation (DNS data). 
Times same as given in figure 1. 
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Fig. 6c: Time t3 



Fig. 6d: Time t4 



Figure 6: Comparison between DNS data and model for 5 ><c (fc). The cut-off wavenumber 
A c = 12. Times of comparison are the same as given in figure 1. 
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